Dynamical analysis of a stochastic delayed SIR epidemic model with vertical transmission and vaccination

In order to describe the dynamic process of epidemic transmission with vertical transmission and vaccination in more detail and to better track the factors that lead to the occurrence of epidemics, we construct a stochastic delayed model with a specific functional response to describe its epidemic dynamics. We first prove the existence and uniqueness of the positive solution of the model. Moreover, we analyze the sufficient conditions for the extinction and persistence of the model. Finally, numerical simulations are presented to illustrate our mathematical findings.


Introduction
A mathematical model has always been an important tool in the study of infectious diseases; there are many works about epidemic such as SIS, SIR, SEIR, and so on. In recent years, research on infectious diseases has been developing, and some good results have been achieved. Elaiw and Agha considered the delayed partial differential equation model to analyze Oncolytic virotherapy based on previous work about ODE [1]. In [2,3], they all applied different models to analyze the SARS-CoV-2 and some useful suggestions were put forward. There are many other papers about epidemic model, such as [4][5][6]. All of these works show that it is an effective method to analyze infectious diseases using the infectious disease model. However, in real life, some infectious diseases may be transmitted vertically from one person to another; that is, the offspring of infected parents may be infected with infectious diseases such as hepatitis and tuberculosis at birth, called vertical transmission [7], so how to effectively prevent and control the spread of infectious diseases has become an important topic in epidemiology, on the current research shows that vaccination has become an important and commonly used strategy to eliminate infectious diseases, it can effectively reduce the infection of infectious diseases [8][9][10]. Based on the SIR epidemic model with vaccination and vertical transmission model proposed by Meng and Chen [11], we establish the following model, and a framework diagram of the disease model is shown in  In the above model, S(t), I(t), and R(t) denote the number of susceptible, infective, recovered individuals at time t, respectively. We suppose that b and b are the birth rate coefficients of the non-infected person (S + R) and infected person; I, d, and d are their corresponding death coefficients, respectively. The infection rate of the disease is β; the susceptible person is an infected as infected person at a bilinear rate of βI(t), and the infection recovery rate is γ . The proportion of the offspring of infectious parents who are susceptible is p; the proportion of the offspring of infectious parents who are infected is q, 0 < p < 1, 0 < q < 1 and p + q = 1. m (0 < m < 1) is the proportion of the successfully vaccinated population to the entire susceptible population.
But in fact, the vaccine validity is usually limited, and the immunized person's immunity can disappear [12][13][14][15][16]. Suppose ω denotes the period of vaccine validity, then the susceptible inoculated at tω will become susceptible again at t. However, due to the existence of natural mortality, the probability of these vaccinated people still alive at t is e -dω [17]. Then, we can get the model as follows: (1.1) Here, we assume b, d, b , and d to be equal; the case can be seen in [18,19]. It can be found from the model that d(S(t)+I(t)+R(t)) dt = 0 and the population has a constant size, which is normalized to unity. Our analysis below is simplified with this assumption. By calculations, the basic reproduction number of model (1.1) is obtained R 1 0 = βb (b+m(1-e -bω ))(pb+γ ) , and we find that when R 1 0 < 1, the model has a disease-free equilibrium point (S 0 , 0, R 0 ), where b+m(1-e -bω ) , when R 1 0 > 1, the model has an endemic equilibrium point (S e , I e , R e ), where S e = bp+γ β , I e = b b+r (1 -1 . To model the disease transmission process, several authors improved the following bilinear incidence rate βSI to get a more suitable infection rate, where β is a positive constant [20]. However, there are many forms of nonlinear incidence rates for a more generalized form. The incidence rate βSI/(1 + α 1 S + α 2 I + α 3 SI) was introduced by Hattaf et al. [21] and used in [22]; it is a general form, which represents mutual interference between S and I, βI and measures the infectivity of the disease when it enters a fully susceptible population. 1/(1 + α 1 S + α 2 I + α 3 SI) measures the inhibition effect from the behavioral change of the susceptible population and the infected population when their number increases or from the crowding effect of the infected individuals, that is, due to the information of the disease. The infected or susceptible population will cause behavioral changes and inhibitory effects; therefore, it is more interesting and valuable than bilinear incidence rate. It has been widely applied in epidemiological studies. When α 1 = α 3 = 0, then we have the saturated incidence rate βSI/(1 + α 2 I) [23], which was used in [24][25][26]; when α 3 = 0, we get the Beddington-DeAngelis functional response βSI/(1 + α 1 S + α 2 I) [27], which was used in [28,29]; when α 3 = α 1 α 2 , we get the Crowley-Martin functional response βSI/(1 + α 1 S + α 2 I + α 1 α 2 SI) [30]. α 1 , α 2 , and α 3 are used to measure the inhibitory effect on infectious diseases when crowding effect or behavioral changes caused by the increase of susceptible individuals occur, and the infection coefficient can be effectively reduced by reasonably selecting appropriate parameters.
In this article, a Crowley-Martin functional response is considered; that is, the incidence rate of disease is modeled by βSI/f (S, I) = βSI/(1 + α 1 S + α 2 I + α 1 α 2 SI), where β is the infection coefficient and α 1 , α 2 ≥ 0 are constants. Thus, we get (1.2) For model (1.2), we can use the results presented by Hattaf et al. [31]. It is easy to get the basic reproduction number of disease that is given by On the other hand, the current environmental fluctuations have a great impact on all aspects of real life, so we will consider the impact of environmental fluctuations on the transmission rate β. Unless otherwise specified, it is assumed here that the random disturbance is a type of white noise, namely β dt → β dt + σ dB(t), where B(t) is a Brownian motion, and σ is intensity. We let B(t) be defined on a complete probability space ( , F , {F t } t≥0 , P) with a filtration {F t } satisfying conditions that are increasing and right continuous while F 0 contains all P-null sets. Then the form of the stochastic model corresponding to the deterministic model (1.2) is as follows For biological significance, the following analysis satisfies the condition S(t) ≥ 0, I(t) ≥ 0, R(t) ≥ 0. Then, noticing the first two stochastic differential equations in system (1.4) do not depend on the function R(t), we can exclude the third one without loss of generality [32,33]. Hence, we will only discuss the following system: (1.5) This paper is organized as follows: in Sect. 2, the global existence, positivity, and boundedness of solutions of our stochastic model (1.5) will be proved. In Sects. 3 and 4, we respectively show sufficient conditions for the extinction and persistence of the disease.
In Sect. 5, some numerical simulations are presented to illustrate our main results. Finally, the paper ends with a brief discussion and conclusion in Sect. 6.

Existence of the positive solution
In this section, we establish the global existence, positivity, and boundedness of solutions of system (1.5). Since S(t) and I(t) in system (1.5) denote population sizes, they should be nonnegative, so for further study, we should firstly give region to prove that system (1.5) has a unique global positive solution. First, we can find that it is clear that , Proof Let (S(θ ), I(0)) ∈ , θ ∈ [-ω, 0) and n 0 > 0 be sufficiently large such that each component of (S(θ ), ]. Define, for each integer n ≥ n 0 , the stopping times , It suffices to prove that P(τ = ∞) = 1, that is P(τ < t) = 0, ∀t > 0, we can see clearly that We only need to show that lim sup n→∞ P(τ n < t) = 0 for this matter, referring to [35,36], we can take a similar function and use some approaches to prove the theorem. Then, we set a C 2 -function U: R 2 + → R + for all (S(t), I(t)) > 0: Applying Itô's formula, for all t ≥ 0 and s ∈ [0, t ∧ τ n ], we obtain dU X(s) Then Taking integral and expectations on both sides of (2.1) and applying Fubini's theorem, we get Using Gronwall's inequality, we have Hence, Since U(X(t ∧ τ n )) > 0 and some component of X(τ n ) is less than or equal to 1 n , we deduce that Thus, lim sup n→∞ P(τ n < t) = 0. This completes the proof.
The following theorem proves that there is a unique globally positive solution to system (1.5) for any initial value X 0 = (S(θ ), Theorem 2.2 For any initial value S(θ ) ≥ 0 and I(0) > 0, ∀θ ∈ [-ω, 0), system (1.5) has a unique positive solution (S(t), I(t)) on t > 0, and the solution will remain in R 2 + with probability one, that is to say (S(t), I(t)) ∈ R 2 + for all t > 0 almost surely.
Proof Since the coefficients of system (1.5) satisfy the local Lipschitz conditions, then for any initial value S(θ ) ≥ 0 for all θ ∈ [-ω, 0) and S(0) > 0, I(0) > 0, there is a unique local solution (S(t), I(t)) on t ∈ [0, τ e ), where τ e represents the explosion time [37]. To verify this solution is global, we only need to show τ e = ∞ a.s. To this end, let k 0 ≥ 1 be sufficiently large such that (S(θ ), I(0)) belongs to the interval [ 1 k 0 , k 0 ]. For each integer k ≥ k 0 , let us define the following stopping time where throughout this paper, we set inf ∅ = ∞ (as usual, ∅ represents the empty set). Obviously, τ k is increasing as is true, then τ e = ∞ a.s. and (S(t), I(t)) ∈ R 2 + a.s. for all t > 0. That is to say, to complete the proof, we only need to show τ ∞ = ∞ a.s. If this assertion is false, then there exists a pair of constants T > 0 and ∈ (0, 1) such that P{τ ∞ ≤ T} > . Thereby, there is an integer k 1 ≥ k 0 such that the nonnegativity of the above function can be seen from u -1 -ln u ≥ 0 for ∀u > 0, let k ≥ k 0 and T > 0 be arbitrary. Applying Itô's formula, for all t ≥ 0, we can get: where , here because f (S, I) ≥ 0, so LV (S, I) ≤ β + 2b + m + pb + γ + σ 2 =: K , where K > 0 is a constant. So, the above can be written as Set k = {τ k ≤ T} for k ≥ k 1 and by virtue of (2.4), we obtain P( k ) ≥ . Note that for every ω ∈ k , the S(τ k , ω) or I(τ k , ω) equals either k or 1 k . Consequently, V (S(τ k , ω), I(τ k , ω)) is no less than either k -1 -ln k or 1 k -1 -ln 1 k = 1 k -1 + ln k. Hence, we can get It follows from (2.6) that where I k is the indicator function of k . Letting k → ∞, then we have ∞ > V (S(0), I(0)) + KT = ∞, which yields the contradiction, we get τ ∞ = ∞. This means that S(t) and I(t) will not explode in a finite time, almost surely. This completes the proof.

Extinction of the disease
In this section, we study the extinction of the disease. Before giving our main result of this section, let us present some lemmas. Proof The proof is similar to that in Liu et al. [38] and hence is omitted.
Proof It follows from Itô's formula that integrating this from 0 to t and dividing by t on both sides, we have where M(t) = taking the limit superior of both sides, we obtain the desired assertion (3.1).
On the other hand, noting that the Crowley-Martin functional response can be written differently as By the law of large number for martingales and for R s 0 < 1, we obtain ln I(t) t ≤ R 2 0 1 - taking the limit superior of both sides, we get the assertion (3.2). We have proved that holds. This completes the proof.

Persistence of the disease
In this section, we study the conditions for the Persistence of the disease. For simplicity, we define X(t) = 1 t t 0 X(s) ds and present the definition of persistence in the mean as follows [37]. where Proof Since (S(t), I(t)) ∈ , from the Crowley-Martin functional response (3.4), we get I.
integrating both sides of (4.3) from 0 to t, there is Thus .
In view of Lemma 3.1, one can easily obtain that by (4.4) and (4.6), we get = 0 a.s., by Lemma 4.1 and R s 0 > 1, we deduce that This is the required inequality (4.1), and from (4.5), we have Therefore, I * > 0.

Simulations
In this section, we will use the Milstein method and the Euler-Maruyama method [40] to illustrate our results, and all the step sizes are 0.1 [41]. We take 50 realizations and use their average to plot such as I(t) = 50 i=1 I i (t)/50, where the I i (t) represents the ith realization. We compare the threshold parameters of the deterministic model and stochastic model to explain the effect of white noise on the system. A typical example of vertical contagious and vaccine-related infectious diseases is hepatitis B. There are many studies on hepatitis B. In this part of numerical simulation, the value of parameters is taken from [42][43][44][45]. For the stochastic model (1.5), we consider the discrete equation: Here ξ k (k = 1, 2, . . .) is the N(0, 1)-distributed independent Gaussian random variables. Now σ (t) is the intensity of white noise and time increment t > 0. For the deterministic system (1.2), we choose the initial value (S(0), I(0)) = (0.5, 0.2) and the parameter values β = 0.4, b = 0.3, p = 0.1, m = 0.9, α 1 = 0.6, α 2 = 0.1, γ = 0.2. We compare the two cases of ω = 1 and ω = 2 when β = 0.4. By simple calculations, we get both of them R 2 0 < 1 and find that the I tend to 0, which means that the disease fades. In Fig. 2(a), (b), we see that the larger the value of the time delay ω, the faster the disease will fade. When β = 0.6 and ω = 1, we get R 2 0 > 1, it shows that the disease becomes endemic, but when we increase the time delay to 2, we get R 2 0 < 1 and find that the disease is extinct from From the above analysis, it is found that the longer immune period, that is, time delay ω, the less likely the disease will break out. For the stochastic system (1.5), we first use the Milstein method [40] to illustrate our results, and we choose the initial value (S(0), I(0)) = (0.7, 0.4) and the parameter values b = 0.3, p = 0.5, m = 0.2, α 1 = 0.6, α 2 = 0.1, γ = 0.2, ω = 1. When β = 0.4 and σ = 0.5, we get R s 0 = 0.5310 < 1 and σ 2 ≤ β(b+α 1 b+m(1-e -bω )) b ; hence, the condition (b) of Theorem 3.1 is satisfied.
When β = 0.4 and σ = 0.9, we get σ 2 > β 2 2(pb+γ ) , the condition (a) of Theorem 3.1 is satisfied. In Fig. 3(a), (b), the I both exponentially decays to zero, which indicates the extinction of the disease. Next, we let parameter β = 0.8 and σ = 0.5 and others are the same as above. In this case, we get R s 0 = 1.2712 > 1, according to Theorem 4.1, the disease is persistent, see Fig. 3(c). As shown in Fig. 3(d), when σ increases to 0.9, the I also exponentially decays to zero, which indicates the extinction of the disease. The above results show that, to a certain extent, stochastic noise has an effect on infectious diseases, properly increasing noise intensity can reduce the spread of infectious diseases.
On the basic, we use Euler-Maruyama method and present the parameter values β = 0.7, b = 0.3, p = 0.5, m = 0.2, α 1 = 0.6, α 2 = 0.1, γ = 0.2. Here we analyze the impact of time delay ω changes on infectious diseases when the noise intensity σ is little, i.e., σ = 0.2. We find that when ω = 1, the disease is persistent and is shown in Fig. 4(a), when ω = 2 the disease is beginning to go extinct, until ω = 3 or ω = 5, the disease has become extinct, and they are shown in Fig. 4(b), (c), (d), respectively. Based on the above analysis, we can get Figure 4 The effect of delay ω on dynamics of the stochastic system (1.5) that the time delay can contribute to the extinction of the disease, the larger the value of the time delay ω, the faster the disease will fade.
Next, we use the Milstein method [40] to analyze the influence of some parameters. We choose the initial value (S(0), I(0)) = (0.5, 0.4) and the parameter values β = 0.8, b = 0.3, m = 0.2, α 1 = 0.6, α 2 = 0.1, γ = 0.2, σ = 0.2, ω = 1. Under these parameter values, we choose that in contrast, p = 0.4, p = 0.6 and p = 0.8, see Fig. 5(a). The higher the p value, the more susceptible, the fewer people vertically infected; therefore, the disease will be quickly controlled. Then, we consider the proportion of successfully vaccinated population m, choose p = 0.5, and other parameters do not change. Taking different values of m, we get that the bigger proportion of the successfully vaccinated population, the less infected, see Fig. 5(b).
Lastly, we briefly describe the immunity level. We consider the form of successful immunization as dV (t) dt = mS(t) -mS(tω)e -bω and analyze how the immunity level changes as the time delay ω and the decay rate b change. We set β = 0.8, p = 0.5, m = 0.3, α 1 = 0.6, α 2 = 0.1, γ = 0.2, σ = 0.1, and the initial value (S(0), I(0)) = (0.5, 0.4). Then, we can find when b is larger, the V ((t)) is higher in Fig. 6(a). The higher the decay rate b, the smaller the probability of death due to disease, and then, the population size of successful vaccination is larger. On the other hand, we set b = 0.3, and other parameters do not change.

Figure 5
The effect of some parameters on dynamics of the stochastic system (1.5) Figure 6 The immunity level Then, in Fig. 6(b), we can find when ω is larger, the V ((t)) is higher. The greater the time delay ω, the less likely it is to die from the disease, and then, the population size of successful vaccination is larger. This shows that the effective period of vaccine is the longer, the immunity level is higher.

Conclusions
In this paper, we have analyzed a stochastic delayed SIR epidemic model with vertical transmission and vaccination, the introduction of stochastic effect and time delay into deterministic models gives us a more realistic way of constructing epidemic model. In addition, we consider a specific functional response incidence rate. In this model, firstly, we have proved the global existence, positivity, and boundedness of the solution. In addition, we have shown that the disease fades when the white noise is large enough such that σ 2 > β 2 2(pb+γ ) . Moreover, when the noise is small, i.e., σ 2 ≤ β(b+α 1 b+m(1-e -bω )) b , the extinction of the disease can be determined by the value of R s 0 ; if R s 0 < 1, the disease fades. The persistence of the disease is determined by R s 0 , i.e., if R s 0 > 1, the disease persists. Finally, we have simulated our theoretical result and have also found that when the white noise is small, the stochastic system is similar to the deterministic system, but when the white noise is large enough, the stochastic system will appear to be a different phenomenon. Large white noise can suppress the spread of disease. And the higher the p value, the higher the proportion of susceptible newborns and the fewer patients with vertical transmission. With the development of modern medical treatment, there will be more medical measures to block the vertical transmission of infectious diseases. In addition, the increase in vaccination rate m and time delay ω have some influence on the development progress of the disease, they can effectively suppress the occurrence of the disease under the right circumstances.
It is very meaningful to study the epidemic model with time delay caused by vaccination, and we can not ignore the influence of vaccination on some infectious diseases. Moreover, by studying the dynamic behavior of stochastic infectious disease system, we can reflect the actual phenomenon more accurately and reveal the influence of stochastic disturbance on infectious disease system, which is of great significance for the scientific prediction of disease development trend and epidemic prevention and control. It can help us offer some useful control strategies to regulate disease dynamics. In future work, we will consider the delayed SIR model with different incidences to build more realistic models and analyze some other characteristics about them.